function loga()
clc;

loga2(8101, 6, 7531);
loga2(29, 2, 18);
loga2(11,4,9);
loga2(101,13,17);
end

function result = loga2(p, a, beta)

calc = p - 1;
res = zeros(1);
i = 2;
j = 1;
k = 1;
% tozstenyezoi
while (k < 100)
    if ( mod(calc,i) == 0)
        res(j) = i;
        calc = calc / i;
        j = j + 1;
    else
        i = i + 1;
    end
    if calc == 1
        break;
    end
    k = k + 1;
end

uniq = unique(res);
nr = size(uniq);
nr = nr(2);

r = -1*ones(4);
c = -1*ones(4);

for j=1:nr
    q = power(uniq(j), res);

    kpos = 1;
    temp = beta;
    for i = 1:q
       
        if (i == 1)
            r(uniq(j),i) = modmod(temp, ((p-1)/uniq(j)), p);
        else
            r(uniq(j),i) = modmod(temp, ((p-1)/(uniq(j)^i)), p);
        end
       
        for k = 0:uniq(j)-1
            if ((r(uniq(j),i)) == modmod(a, mod(k*(p-1)/uniq(j),p), p))
                c(uniq(j),kpos) = k;
                kpos = kpos + 1;
                break;
            end
        end
        c_temp = c(uniq(j), kpos-1) * (uniq(j)^(i-1));
        temp = mod(temp*modmod(a, mod(p-c_temp-1,p), p), p);
    end
end

% r
% c

M = zeros(1,nr);
X = zeros(1,nr);
y = zeros(1,nr);
x_nr = 1;
for i = 1:nr
    q = power(uniq(i), res);
    for j = 1:q
        X(x_nr) = X(x_nr) + c(uniq(i),j) * uniq(i)^(j-1);
    end
   
    powers = uniq(i)^q;
    M(x_nr) = (p - 1) / powers;
   
    y(x_nr) = modmod(M(x_nr), p, powers);
   
    x_nr = x_nr + 1;
end

X
M
y

modmod(4,8,11)

result = 0;
for i = 1:nr
    result = result + X(i) * M(i) * y(i);
end

result = mod(result,p-1);
fprintf('p = %d a = %d b = %d\n',p,a,beta);
fprintf('%d^%d mod %d = %d\n',a,result,p,beta);
modmod(a,result,p)

end

function c = modmod(number, power, p)
    c = 1;
    for i = 1:power
        c = mod(c*number, p);
    end
end

function c = power( number, res )
    nr = size(res);
    nr = nr(2);
    count = 0;
    for i = 1:nr
        if res(i) == number
            count = count + 1;
        end
    end
    c = count;
end